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Abstract This is the first survey of highly excited eigenstates of a chaotic 3D billiard. We 
introduce a strongly chaotic 3D billiard with a smooth boundary and we manage to calculate 
accurate eigenstates with sequential number (of a 48-fold desymmetrized billiard) about 45,000. 
Besides the brute-force calculation of 3D wavefunctions we propose and illustrate another two 
representations of eigenstates of quantum 3D billiards: (i) normal derivative of a wavefunction 
over the boundary surface, and (ii) ray — angular momentum representation. The majority 
of eigenstates is found to be more or less uniformly extended over the entire energy surface, 
as expected, but there is also a fraction of strongly localized — scarred eigenstates which are 
localized either (i) on to classical periodic orbits or (ii) on to planes which carry (2+2)-dim 
classically invariant manifolds, although the classical dynamics is strongly chaotic and non- 
diffusive. 
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The structure of individual quantum eigenstates of classically chaotic closed systems is one of 
the most important and difficult questions of quantum chaos. Understanding of localization 
properties of highly excited eigenstates can explain almost any statistical property of a quantum 
system, such as e.g. statistics of energy levels (see e.g. ]l]]) or statistical distribution of transition 
amplitudes . For classically fully chaotic — ergodic billiards we have a theorem ||] which states 
that asymptotically, in the semiclassical limit h — ► 0, probability density of eigenstates should 
become uniformly delocalized — extended over the billiard domain, i.e. phase space distributions 
of eigenstates should approach a uniform distribution over the energy surface. The fluctuations of 
an eigenfuction \P(r) are also expected to obey a universal law for classically fully chaotic systems 
Q. For systems having a time-reversal symmetry, ^(r) is expected to behave like a real Gaussian 
random variable with a uniform probability distribution ctP/d^ = l/(\/2vro") exp(— ty 2 /2a 2 ) 
where inverse variance 1/cr 2 is constant and equal to the volume of the billiard. 

For finite values of Planck's constant fi there may be considerable deviations from this 
semiclassical limit: (i) There may be classical adiabatic invariants, which change slowly and dif- 
fusively in the course of time, or partial barriers in phase space — cantori, so that classical orbits 
spend a long time to uniformly fill entire energy surface. If this time is longer than Heisenberg 
(break) time tbreak = frd(E), where d(E) is the density of states, then one should expect strong 
localization of quantum eigenstates on such approximately disconnected components of energy 
surface Q. (ii) But even if classical dynamics quickly explores the entire energy surface w.r.t. 
tbreak one can get strong enhancements of probability density of eigenstates on the least unstable 
classical periodic orbits — the so-called scars ||. 

Semiclassical theory has been developed, which explains the effect of classical periodic orbits 
onto the composition of eigenstates within an energy interval much larger than the mean level 
spacing ||, but there is no adequate theory to predict the structure of individual eigenstates. 
One can use a heuristic argument due to Heller Q , saying that an eigenstate is scarred by a given 
classical periodic orbit if a Gaussian wavepacket, which is launched along the orbit, interferes 
constructively with itself after one period. But in chaotic systems there are many periodic orbits 
and it is typically impossible to tell afortiori which orbits will scar most intensively. (The rule 
that these are the least unstable ones holds only on average.) 

In the literature there exists a vast amount of numerical evidence on scarring in chaotic 
systems with two freedoms (see e.g. JrfJ , |§ , ) , while the structure of comparably high- lying 
eigenstates of three and higher dimensional chaotic systems have numerically not yet been stud- 
ied, so it is not known whether classical periodic orbits can scar eigenstates in higher dimensions 
or not. Heuristic geometrical argument might suggest that scarring by individual periodic or- 
bit is less effective in 3D than in 2D since more phase space is available for a wavepacket to 
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disperse before it can interfere with itself. But scarring by families of similar periodic orbits 
may be effective in higher dimensions. Actually, as approaching the limit Ti — > 0, quantum 
invariant states — eigenstates should approach classical invariant states which are spanned by 
characteristic functions on the classically invariant (sub)manifolds in phase space; but only the 
microcanonical distribution is stable for long times in a fully chaotic system, so it is the only 
one which survives the semiclassical limit. Periodic orbits are only special cases of classically in- 
variant submanifolds in the phase space of a dynamical system ((1+1)D invariant submanifolds 
when they are parametrized by the energy E). In a 3D system there may exist also (2+2)D 
invariant submanifolds in 6D phase space which can have 2D (or 3D) projection onto 3D con- 
figuration space supporting a 2D (or 3D) scar. This type of localization is a genuine 3D effect 
which is yet to be numerically confirmed. 

The purpose of this paper is to give the first survey of the high-lying eigenfunctions of a strongly 
chaotic 3D system. We have defined a family of generic 3D billiards with a smooth C°° bound- 
ary (see also [pi]] ). Since no generic system is known to be rigorously ergodic we have defined a 
two-parameter family of 3D billiards whose shapes are given by simplest smooth deformations 
of a sphere: The radial distance re(fi) from the origin to the boundary as a function of the 
direction n,n 2 = 1, is 



and contains the two lowest order terms which preserve the cubic symmetry (the first two fully 



exploration of a parameter space (a, b) we have decided to chose: a = —l/5,b = —12/5. For 
these values of the parameters the classical billiard is strongly chaotic: nondiffusive, without 
partial barriers in phase space, and with large average maximal liapunov exponent {\rnax 

) = 0.54 

meaning that roughly after five bounces we loose one digit of information of a classical orbit. 
There is a tiny regular component of phase space whose relative volume has been estimated to 
be p\ ~ 10 -3 . The billiard domain is marginally convex and there are few isolated neutrally 
stable — parabolic periodic orbits which touch the boundary at the points of zero curvature 
radii. 

The billiard is invariant under 48- fold cubic symmetry group Oh- So we desymmetrize it: 
In the following we consider a 1/48 of an original billiard ([l]) with x > y > z > 0, i.e. we 
consider a 3D billiard in a domain, bounded by the boundary surface B = {r;rs(r/r) = r}, 
and the symmetry planes z = 0, x = y, and y = z. The symmetry planes in phase space, 
fcz = {z = 0,p z = 0}, Q xy = {x = y,p x = p y }, and Q yz = {y = z,p y = p z }, and the boundary 
surface with the tangential momentum S = {r^(r/r) = r,p- V(r^(r/r) — r) = 0} are invariant 






cubic harmonics after Von Lage and Bethe |10| ). After a careful numerical 
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(2+2)D submanifolds with respect to the classical dynamics which have 2D projections onto 3D 
configuration space. When b = 12a, as is the case with our choice, the billiard (j^) possesses 
another less trivial example of 4D invariant submanifold, namely T = {x = y + z,p x = p y +p z }, 
with a 2D configurational projection, the plane x = y + z. 

Quantum eigenfunctions of our billiard (|l]) carry ireducible representations of a cubic group. 
We have decided to study only the singlet eigenstates belonging to fully symmetric ID irrep. 
This corresponds to a study of a desymmetrized billiard, Helmholz equation (V 2 + k 2 )fy ^(r) = 0, 
with Von Neuman boundary conditions on the symmetry planes z = 0, x = y and y = z, and 
Dirichlet boundary conditions on the boundary surface B. Quantum eigenfunctions $(r) have 
been calculated using 3D generalization of scaling method, proposed recently by Vergini and 
Saraceno |]l2[ . For billiards whose shape is geometrically not very far from a sphere, we propose 

to use spherical waves 4>kim(r) = ji{kr)Yi m (r/r) insetead of plane waves exp(i/c • r) as the basis 
of scaling functions which already solve the Helmholz equation but do not satisfy the boundary 



conditions. We define the efficiency of the basis of spherical scaling functions [11] as rj = V/V sp h 
where V is the volume of the billiard, and V sp h is the volume of the smallest sphere enscribed 
to the billiard which has radius r max = max{rB(r/r)}. For a desymmetrized billiard one should 
use only the linear combinations of spherical waves which already satisfy Von Neuman boundary 
conditions on the symmetry planes. These are the products of spherical Bessel functions and the 
fully symmetric cubic harmonics (of type (a) after Von Lage and Bethe ji~0|])) — cubic waves, 
labeled by a, 

<Pka(r) = ji a (kr) c *j{Yi a Aj{r/r) + Y la _ 4j (f/r)). 

3=0 

The coefficients c a j are determined by the Gram-Schmidt orthogonalization of the columns of 
the projector onto the fully symmetric irrep of the cubic group Oh,P = (1/48) J2g£O h in 

the basis of spherical harmonics Y\ a ^j with even angular momentum l a = 2n, Pj"/ = (1 + 

2D( 2n )(0,7r/2,0) 4ji4i /)/3 (using the standard notation 0). 

For an accurate quantization at fixed wavenumber k we need cubic waves a with angular 
momentum l a up to 

Irnax(k) ■ — kT max + evanescent- 

The number of included evanescent angular momenta Al evanescen t is found to be proportional to 
the number of accurate digits of the results, and being much smaller than kr max . For any even 
I, I < Imax there are |_Z/12j + 1 — 5 r ^2 fully symmetric cubic waves, where I = 12|_//12j + r, < 
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r < 12, which sum up to the total number 

N CW (k) = l max (fe) 2 /48 + (fc)/4 + 0(l) 

of fully symmetric cubic waves </>fc a (r) which should accurately capture the eigenfunction ^(r) 
at the wavenumber k 

^k(r) =J2^ka{r)- (2) 

a 

One should find such k and coefficients ip a that ^f(r) satisfy the boundary condition on the 
boundary surface , = 0. Following |12| we minimize a special boundary norm 



f(k) = I #||^(r1| 2 , (3) 



v ■ r 



B 



where d 2 S is a boundary surface element and V = V(rs(r/r) — r)/| . . . | is a unit vector normal 
to the boundary B, by solving the following generalized eigenvalue problem 

~jj: Faa ' ( fc o)Vv = A ^ Oo)Vv , (4) 

a' a' 

where 

f d 2 S 

B 

is a positive definite matrix Ncw(ko) x -^cw(^o)- Every matrix element is an integral over 
the surface of the a sphere (or over 1/48 of it), with integration measure written in spherical 
coordinates (9, </>) 

d 2 S 
v ■ r 



d cos 9d(f) J r\ + (<%r B ) 2 + (c^re/sim 



where tb = re(sin # cos 0, sin sin 0, cos#), which is accurately numerically evaluated using a 
product grid of a Gaussian quadrature for cos# integration Jq d cos 9 with l ma x/^ (positive) 
nodes and a uniform grid for (p integration Jq^ 4 deft with l m ax/^ nodes. [] The eigenvalues A*™) of 
largest modulus are good numerical estimates of billiard's eigenvalues k^ ~ ko — 2/\( n \ Due 
to the scaling symmetry, </>fc iQ (r) = <pk ,a(kr/ko), the corresponding eigenvectors ip^ yield very 

accurate eigenfunctions $W(r), via expansion (|2|), with error being of the order 0((k^ — feo) 4 )- 

2 We have also implemented an optimal cubature formula due to Lebedev for accurate numerical integration 
over the sphere, which reduces the size of a grid for integration over the desymmetrized boundary surface B, 
w.r.t. product grid, by a factor of 9/2, and the total numerical labour by 20 — 40%. Unfortunately we found 
that numerical computation of the grid points and weights of such an optimal grid for large l ma x turns out to be 
extremely difficult numerically, more than the problem itself. 
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Accurate eigenvalues have been obtained by minimizing the boundary norm f(k) (as expanded 
in a Taylor series around ko up to 8th order) for fixed coefficients ifta . The number of converged 
eigenvalues k^ or eigenvectors in a single diagonalization of (0) depends on the required 
accuracy, empirically we find |llj] that it is proportional to k 2 . 

We have chosen ko = 360.0 and calculated a stretch of 64 consecutive eigenstates with eigenval- 
ues k^ from a narrow window centered around ko, corresponding to sequential quantum number 
of a desymmetrized billiard M(ko) = 45, 103 ± 2 (using 3D Weyl formula [0]). Computation 
has been performed with levanescent = 40 (and l ma x = 364) and the boundary norm ([|) of nor- 
malized wavefunctions ty( n ) has ranged between 5 • 10~ 9 and 5 • 10~ 8 . However, it is much easier 
to calculate high-lying eigenfunctions of a 3D billiard, than to analyze and faithfully represent 
them without loosing much of information (since we can plot functions of at most two variables). 
For the presentation we have picked a sample of 12 eigenstates out of 64: 11 consecutive states 
and a specially chosen state with the criterion of being most intensively localized. 

The first representation: normal derivative on the boundary surface V ■ V^"^. Although 
this representation formally contains all the information about a given eigenstate, it explicitly 
displays only a very limited portion of configuration space and it is not very suitable to detect 
localization which can take place strictly inside the billiard region. Nevertheless, it is the easiest 

to be calculated numerically, v ■ V\l/( n ) = Sa^cT V<p k ( n ) a , one only has to calculate the 

normal derivatives of the cubic waves on B. In figure 1 we show the nodal lines of V ■ V^ 71 ) 
on B and the contour plots of its magnitude. Actually, for easier planar presentation we map 
the desymmetrized boundary B, which has a shape of a deformed spherical triangle with angles 
7r/4,7r/2,7r/3, at the points A 1 = (1, 0, 0)r B (l, 0, 0), A 2 = 2~ 1 / 2 (1, 1, 0)r B (2- 1/2 , 2" 1 / 2 , 0), A 3 = 
3~V 2 (1, i ; l)r j e(3^ 1 / 2 , 3 -1 / 2 , 3 -1 / 2 ), where, respectively, 4-fold, 2-fold, and 3-fold symmetry axis 
penetrates £>, onto the half-square {((ft, x),0 <%<(/>< 7r/4} : 

r(4>,X) = n((j),x)rB{n{4>,x)), 

n(cf>,x) = ((l-z(x) 2 )" 1/2 cos^,(l-z(x) 2 )- 1/2 sin0,z( X )), 
z( X ) = (l + ( S m X )- 2 )- 1/2 
Note that this mapping is not far from being area and angle preserving. So we plot 

B^{4>, X ) = v-V^\ F[(t> , x y 

Figure 1 appears similar to the pictures of eigenfunctions of 2D chaotic billiards. On diagrams 
we indicate the projections of known (2+2)D invariant submanifolds: the symmetry planes are 
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mapped on the boundary of the triangle 1Z Z — > \ = 0, Qxy </> = 7r/4, Q yz — > % = 0, while the 
invariant plane T is mapped on the curve sin % = cos <^ — sin (p. Most of states appear more or less 
random (the statistical distribution of normal derivatives over B has been analyzed and found 
to be very close to a Gaussian for most of states). However the state 12 is strongly enhanced 
in the vicinity of point A\, as well as the states 7,8, 10, while the states 7, 10 are significantly 
enhanced in the vicinity of the projected manifold T. 

Uniform representation of 3D eigenfunctions without loosing essential information: ray - 
angular momentum representation. In the momentum — angular momentum representation, 
quantum states are labeled by eigenvalues of the three commuting observables, p 2 , (r x p) 2 , 
and (p x r) z , namely h 2 k 2 , h 2 l(l + 1), and hm, 

r B (n) 

(klm\^) = Jd 2 nY[ m (n) J drr 2 ji{kr)^(r). 
4tt 

Note that eigenstates in the momentum — angular momentum representation, normalized as 
J2im Io° dkk 2 \ (klm\fy( n ^)\ 2 , are strongly localized around k^ in the momentum coordinate k , 
which is classically a constant of the motion. The structure of eigenfunctions going perpendic- 
ularly to the energy shell k ~ k^ is somehow trivial ({klm\^f) ~ 0, if \k — k^\r max 3> 1) and 
certainly less interesting than the structure inside the shell. Therefore we define a ray — angular 
momentum representation of 3D wavefunctions 

00 

iW*] = / dkk 2 \(klm\*)\ 2 , (5) 


which gives the total probability for a particle with arbitrary momentum hk moving along a ray 
with fixed angular momentum (l,m are the natural coordinates on the spherical energy shell 
in the momentum space). Numerical calculations with formula (j|) are very tedious since they 
require calculation of the full momentum — angular momentum representation (klm\^) first. 
Using well known relations for spherical harmonics and Bessel functions one can reformulate eq. 
(|5|) giving more practical expression 




For the eigenstates \I/( n ) of our billiard we have Rj^ = unless I = 2n, m = 4j and = R\ \l m 
(due to the symmetry), and the integration over r can be simplified up to the radius of the largest 
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inscribed sphere, r min = min{r #(f / V)} , 

\ 2 r 

" drr 2 jf(k^r) 




r 



in a -V 



+ — drr 2 
2 J 





- 2 

dMO € {r B (n) - r)Yi Aj {n)^ n \rn) 



The numerical integrals over r have been performed using another (sufficiently dense) Gaus- 
sian quadrature, and smooth (Fermi) approximation for the step function e (x) = 1/(1 + 

exp(— x/e)), where e < 2ir/k^ n \ For chaotic (microcanonical) eigenfunctions, are ex- 

pected to fluctuate around smooth microcanonical classical angular momentum distribution 

^im = {fiftl — \r x p\)S(hm — (fx p)z)) on the classically allowed region of lm where R& is 
nonzero. Strong systematic deviations indicate localization or scarring. However, the statistical 

fluctuations of pi m = Ri^l/Ri m are veT ^ l ar § e - We have numerical evidence for the stationarity 
and Lorentzian tails of the statistical distribution of pi m , dV/dp(p S> 1) = 0(p~ 2 ), which has 
been studied numerically for the full stretch of 64 chaotic eigenfunctions. 

In Figure 2 we show the sample of 12 states in the ray — angular momentum representation, 
which is found to be very sensitive to localization, since the coordinates I and l z = m are 
the 'most adiabatic' observables, (they are the canonical actions for the spherical billiard). 
Because of strong fluctuations we show (above the diagonals I = m with I and m coordinates 
interchanged) also the smoothed ray-angular momentum representation 

^5 = ^2 E 4"i «p[-(0 - l'f + (m - m') 2 )/(2A 2 )], 

I'm' 

with A = 3. Although the system is classically strongly chaotic, having almost GOE spectral 
statistics [|ll]], we see lots of significantly localized eigenstates. Only the states 1, 2, 3, and 9 are 
really uniformly extended, the states 7 and 10 are localized on the (2+2)D invariant manifold 
T, which is mapped on the line m = i/v3, the states 7, 8 and 10 are localized on the symmetric 
(2+2)D invariant manifolds 1Z x ^ y ^ z , which are mapped on the lines m = and m = I. This 
representation is convenient also because classical periodic orbits are sets of points since each 
line segment of an orbit has fixed angular momentum. The state 12 is a strongly enhanced in the 
neighborhood of an unstable periodic orbit of desymmetrized billiard of length / = 10.209 with 
a pair of unimodular, exp(±il.558), and a pair of real, 4.000, 0.250, eigenvalues of a monodromy 
matrix. 
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Brute-force approach: 3D eigenf unctions in configuration space. Finally, we have calculated 
full 3D wavefunctions for the interesting representative eigenstates of our billiard using expan- 
sion (^) on a dense mesh of points in 3D configuration space. At every point in space we have 
then smoothed the probability density l^l 2 by averaging over a ball of radius equal to one de 
Broglie wavelength 2-Tr/k. The regions of high smoothed probability density (being above certain 
threshold) have been plotted in 3D for states 1,12, and 8, in figures 3,4, and 5, respectively. The 
probability distributions of fluctuation amplitudes ^(r) have been calculated and compared to a 
Gaussian The state 1 (figure 3) is shown as an example of a typical chaotic eigenstate which 
is uniformly extended over the billiard domain with a Gaussian amplitude distribution dV /d^. 
The state 12 (figure 4) is shown as an example of a ID scar localized in a neighborhood of an 
unstable periodic orbit. The state 8 (figure 5) is shown as an example of a 2D scar localized in 
a neighborhood of an (2+2)D invariant chaotic submanifold 1Z Z . (Note that the states 7 and 10 
have a similar structure but they are associated with two invariant (2+2)D chaotic submanifolds 
simultaneously, 1Z Z and T). Also the amplitude distribution dV/d^f for the localized states 8 
and 12 is different from a Gaussian: there is a stronger peak for small amplitudes and slowly 
decaying tails for large amplitudes. 

In a conclusion we wish to emphasize the twofold goal of this paper: (i) the survey of very-high 
lying eigenstates of a strongly chaotic and generic 3D billiard which has been made possible 
by (ii) an application of the scaling method for quantization of billiards (which works for the 



Helmholz equation in arbitrary dimension), proposed recently [Q [. In a preceding paper [11] we 
found small but significant deviations of the level statistics from the expected ultimate asymp- 
totics (GOE) which have been explained in terms of localization of eigenfunctions. Indeed, 
besides the majority of extended states with Gaussian random amplitude, we found examles 
of ID scars which are strongly enhanced in the vicinity of short unstable classical periodic or- 
bits, and examples of 2D scars which are strongly enhanced in the vicinity of the projections 
of (2+2)D invariant chaotic submanifolds on to configuration space. This is the first numerical 
study of the structure of eigenstates of higher dimensional (D > 2) chaotic systems in the semi- 
classical regime (i.e. very large sequential quantum number or small effective Planck's constant). 
Although the phase space is dynamically more 'connected' than in 2D systems, and the classical 
dynamics is really uniformly ergodic on a scale of tb rea k, we have found very rich behavior rang- 
ing from Gaussian pseudo-random extended states to various types of localized states which are 
associated with classically invariant structures. This is the result which should motivate further 
theoretical and numerical study of quantum chaos in higher dimensional systems. 
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Figure captions 



Figure 1: We show the normal derivative of the wavefunctions for the sample of 12 eigenstates 
on the mapped boundary surface, the triangle 0<x<</>< 7r /4 (see text). Below the diagonals 
(abscissa <fi, ordinate \) we show the magnitude of the normal derivative (using 10 levels of 
greyness which change when the square of the function changes for a factor 2.3), and above the 
diagonals (abscissa x, ordinate <p) we show the nodal lines of the normal derivative. The full 
curves going through the middle of triangles are the images of (2+2)D invariant submanifold 
T, whereas the images of symmetric (2+2)D invarinat submanifolds are the boundary lattices 
X = 0, X = 4>, and <p = tt/4. 

Figure 2: Same as in figure 2 but for the ray — angular momentum representation of eigenstates 

Rfol bellow the diagonals (abscissa < / < k^r ma x, ordinate < m < m) and their smoothed 

counterparts above the diagonals (abscissa m, ordinate /). The lines m = //a/3 are the 
images of the relevant (2+2)D invariant manifold T and the image of the relevant periodic orbit 
for the state 12 is indicated with small circles. (See text.) 

Figure 3: The smoothed 3D probability density of a chaotic eigenfunction (state 1 on figs. 1 
and 2). The desymmetrized billiard (we plot the curves along its six edges connecting the four 
corners (below), A\ (left), A2 (middle), and A3 (right) as small circles) is sliced with horizontal 
planes v = (y — z)/V2 = nAv, for At; = 0.005 and n = 0, 1 ... 90, into triangular sections. Going 
from the bottom n = to the top n = 90, the regions where the smoothed probability density 
exceeds a threshold, ^\ hr = 1A/V, are plotted and filled with the level of greyness which is 
proportional to the vertical coordinate v (dark is bottom and bright is top — see grey scale on 
the right). Hopefully, one can thus get an impression of the 3D structure of eigenfunctions. The 
inset on the lower right corner of the picture shows the histogram of the amplitude distribution 
dV/d^ which is in this case in agreement with a Gaussian (dotted curve; vertical dotted lines 
indicate the positions of ±1 standard deviation). 

Figure 4: The same as in figure 3 but for a strongly localized — scarred eigenfunction (state 
12 on figs. 1 and 2). The state is strongly enhanced in the vicinity of an unstable periodic orbit 
(see figure 2 and text), so we have picked a higher threshold Vl/f^ = 1.55/V, and the amplitude 
distribution largely deviates from a Gaussian. 
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Figure 5: The same as in figure 3 but for a strongly localized eigenfunction (state 8 on figs. 1 
and 2) and the same threshold ^1 hr = 1.55/V as in fig. 4. The state is strongly enhanced on the 
projection of a (2+2)D invariant chaotic submanifold 1Z Z , z = 0, so it is a beautifull example of 
a 2D scar. Correspondingly, the amplitude distribution deviates from a Gaussian. 
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